Method of processing seismic data

ABSTRACT

A method of processing seismic data using a seismic energy propagation model of the subsurface is disclosed. The method assigns seismic source, seismic receiver, and reflection point locations to the propagation model; identifies a plurality of alternative raypaths consistent with the propagation model that originate at said seismic source location, reflect at said reflection point location, and terminate at said seismic receiver location; selects a raypath from the plurality of alternative raypaths having a shortest ray length, and utilizes the selected raypath in subsequent seismic processing.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to a method of processing seismic data.

2. Description of Related Art

Seismic data is collected in order to analyse the sub-surface of theearth, in particular for hydrocarbon exploration. Seismic data foranalysing sub-surface structures may be collected on land or, overwater, using sea-going vessels. In order to obtain the data, a seismicsource which may comprise explosives (on land) or an impulse ofcompressed air or airguns (at sea) is provided. The seismic data signalsreflected by the various layers beneath the surface of the earth areknown as traces and are sensed by a large number, typically hundreds, ofsensors such as geophones on land and hydrophones at sea. The reflectedsignals are recorded and the results are analysed to derive anindication of the layer formations beneath the sub-surface. Suchindications may then be used to assess the likelihood of hydrocarbondeposits.

The analysis of the results to derive an indication of layer formations,however, is not straightforward. Particularly where the materials of thesub-surface of the earth vary laterally, there may be more than onesignal path between the seismic source and a point within thesub-surface which reflects the signal. Typically, the same will be trueof the return path between the reflecting point and a respective seismicsensor, such as a geophone or a hydrophone. If a case of three differentpaths in each direction is considered, there will be nine differentround-trip routes by which a signal can travel from the seismic sourceto the seismic sensor from a single reflection point. Cost-effectiveanalysis using all of these possible paths is impossible, so some meansof simplifying the processing is required.

In “Green's Functions for 3D Pre-stack Depth Migration” published in theEAGE 57 Conference and Technical Exhibition in Glasgow, Scotland on 29May to 2 Jun. 1995, the following two prior art techniques for reducingthis complexity are discussed. It will be understood that the signalsresulting from a single reflection point will generally arrive at thegeophone at different times and with different amplitudes dependent uponthe distance of the path travelled and the sound-propagatingcharacteristics of the subsurface layers through which the sound waveshave passed. Consequently, there are a number of “ray paths” through thesub-surface of the earth which relate to a signal reflected by a singlereflection point. One proposed solution to the complexity of thenumerous ray paths is to select a so-called first arrival signal. Thiswill be the arrival signal (or “arrival”) corresponding with the fastestpropagating seismic signal. However, one drawback of this technique isthat the first arrival is rarely the strongest signal and often containstoo little energy to provide reliable and accurate analysis. However,the methods for calculating the first-arriving travel time tend to becheaper and simpler than other methods.

Some of these methods are commonly (and confusingly) referred to as“shortest path” methods, although they actually compute a shortesttravel-time path, rather than a shortest physical ray length path.Articles describing this type of method may also be found in Geophysics,Volume 56, No. 1, January 1991, T. J. Moser, “Shortest path calculationof seismic rays”, pages 59–67; Geophysics, Volume 58, No. 7, July 1993,Robert Fischer et al., “Shortest path ray tracing with sparse graphs”,pages 987–996; and Geophysics, Volume 59, No. 7, July 1994, T. J. Moser,“Migration using the shortest-path method”, pages 1110–1120. Referencesto the process of computing the shortest travel-time raypath (ratherthan the shortest ray length raypath), can be found in these articles onpage 59, abstract, line 4; page 987, column 2, line 20; and page 1111,column 2, line 37, respectively.

Another prior art technique is to select the arrival having the maximumamplitude. However, selection of this arrival is not necessarilystraightforward because the model of the sub-surface will generally onlybe approximate. The maximum amplitude arrival will only provide the bestsingle arrival as long as the estimates of amplitude are correct.Another difficulty with the maximum amplitude arrival is that the choiceof arrival can switch rapidly back and forth between branches. However,improvements using the maximum-amplitude arrival over use of a firstarrival have been observed in the prior art reference first identifiedabove.

It is an object of the present invention to provide a method ofprocessing seismic data which ameliorates the disadvantages of theseprior art techniques.

SUMMARY OF THE INVENTION

According to the present invention, there is provided a method ofprocessing seismic data using a seismic energy propagation model of thesubsurface, including: assigning seismic source, seismic receiver, andreflection point locations to the propagation model; identifyingalternative raypaths that originate at the seismic source location,reflect at the reflection point location, and terminate at the seismicreceiver location; selecting a raypath having a shortest ray length; andutilizing the selected raypath in subsequent seismic processing.

It has been appreciated that applying the shortest ray-path lengthcriterion will virtually never result in a low amplitude signal. Whilethe physical length of the ray will vary in response to changes in thevelocity model, the physical ray length is less susceptible to suchchanges than is the estimate of the amplitude. The reason for this isthat amplitude is related to the curvature of the rays. Thus smallerrors in the ray paths can produce large errors in the amplitudes. Incontrast, the ray-length is an integral quantity (the integral of arclength along the ray) so it is relatively insensitive to smallperturbations in the ray path. In other words, the choice of aparticular ray using the technique of the present invention is lesslikely to switch rapidly than the maximum amplitude technique. Thisprovides a much more reliable assumption upon which to base subsequentprocessing than either of the prior art techniques.

Further preferred features of the present invention are set out in theattached dependent claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention will now be described, by way of example, withreference to the accompanying drawings, in which:

FIG. 1 shows a seismic source and a plurality of seismic receiversarranged above a number of subsurface layers in the earth;

FIG. 2 shows a single source location and a single image location withseveral ray paths between them;

FIG. 3 shows a two dimensional seismic image of distance against depthfor a particular seismic model using the seismic signal having the firstarrival (i.e. the shortest travel time path);

FIG. 4 shows a seismic image of distance against depth for the seismicmodel using the seismic signal having the maximum amplitude arrival; and

FIG. 5 shows a seismic image of distance against depth for the seismicmodel using the arrival associated with the raypath having the shortestraypath length.

DETAILED DESCRIPTION OF THE INVENTION

In FIG. 1, a seismic source S, for example an explosive, is located onthe surface of the earth together with a plurality of geophones R1 toR5. Typically, there will be hundreds of geophones arranged in a twodimensional or a three dimensional array over the surface of the earth.For simplicity, only five geophones R1 to R5 are shown. Seismic signalpaths (“ray paths”) are shown in broken lines between the source S andthe geophones R1 to R5 via subsurface layers L1, L2, and L3. The raypaths are shown reflecting from a number of horizons H1, H2, and H3.Only some of the ray paths are shown for clarity.

Because the seismic signals will travel at different speeds in thelayers (typically the speed of propagation will increase with anincrease in depth) it is not straightforward to locate the horizons H1,H2, and H3. The data from each horizon is generally separated and aseismic energy propagation model (i.e. an estimate of the subsurfaceacoustic velocities) of the subsurface being analysed is applied toderive an estimate of the actual reflecting (or “imaging”) point on ahorizon beneath the surface. The process of prestack depth migrationinvolves simultaneously adding together multiple samples to increase thesignal to noise ratio of the data (similar to a conventional “stacking”procedure), moving seismic events to compensate for the offset distancebetween the source and the receiver (similar to a conventional “normalmove-out correction” procedure), and moving seismic events to compensatefor inclined seismic reflectors (similar to a conventional “migration”procedure). Prestack imaging of seismic data is usually implementedusing an integral formulation in the time space domain.

While the phrase “reflection point” is used throughout this application,this phrase may be more fully be understood as being an “image point”,i.e. a subsurface location illuminated by the raypath in question.Similarly, the raypath from the source to the image point and from theimage point to the receiver will have an abrupt change of directions atthe image point. While this change of direction is referred to as beinga “reflection” throughout this application, it may also be thought of asbeing “scattered” at the image point. The seismic energy propagationmodel may not incorporate any assumptions regarding reflector dip(inclination) angles in the vicinity of the image point and the“downgoing” ray from the source to the image point is not necessarilyrequired to have a seismic reflector incidence angle that is equal andopposite to the seismic reflector incidence angle of the “upgoing” rayfrom the image point to the receiver.

The image at any point is computed from an integral over a surface inthe prestack data. The integral can be written in the generic form of:Image(x _(i))=∫∫WData(x _(s) ,x _(r) ,t _(s) +t _(r))dSdRwhere x_(i), x_(s), and x_(r) are the image, source and receiverlocations, W (x_(i), x_(s), x_(r)) is a (possibly complex and frequencydependent) weight that is a function of the source, receiver and imagelocation, and t_(s)(x_(i), x_(s)) and t_(r) (x_(i), x_(r)) are thetraveltimes from the source and receiver to the image pointrespectively. The traveltimes and weight functions depend on a modelthat is an estimate of the subsurface properties. The integral isevaluated over all source and receiver co-ordinates. The traveltimesdefine a trajectory in the data over which the integral takes place.

The traveltimes may be calculated by many methods (e.g. finitedifference methods, gridded travel time approximation methods,ray-tracing). Most of these methods solve a high frequency approximationto the wave equation which decouples the solution into two parts, firstsolving the eikonal equation for the traveltimes and then solving thetransport equation for the amplitudes. The eikonal equation is anequation that results from asymptotic expansion of the wave equation. Itis a non-linear differential equation that is satisfied by thetraveltimes. All of the methods either explicitly or implicitlycalculate the ray-path which is the path that the energy travels alongbetween the source or receiver and the image point.

When the present inventive method is used, it is preferable to use aray-tracing method. Ray-tracing methods often inherently calculatetravel distances as they select alternative ray-paths or provide outputsthat allow the travel distances of the various ray-paths to be easilycalculated.

However, in a complex subsurface model there may be several ray pathsthat connect the source to the receiver and this means that thetraveltime functions will be multivalued. FIG. 2 shows one example of acomplex model that gives rise to multivalued traveltimes. A seismicsource S is being used to analyse an image I at a point beneath thesurface 10 of the earth. However, there are two salt bodies 12, 14arranged on either side of the direct path 18 between the source S andimage I. The salt bodies 12, 14 refract the seismic rays and provide twofurther paths 16, 20 between the source S and the image I. The two saltbodies thus produce refracted ray paths from the source to the imagepoint in addition to the direct ray path that does not travel throughthe salt.

When the traveltime is multivalued, the correct way to evaluate theintegral is to sum over all branches of the traveltime functions.However, if there are three paths from the source S to the image I andthree from the image to the receiver, then there are 9-branches overwhich to sum. 3-D prestack depth migration is already expensive. A9-fold increase in computational complexity will severely lengthen theanalysis and hence increase the cost greatly. As the number of pathsincrease, the complexity increases accordingly.

If we do not use the results from all of the branches, we must make achoice of which branch or branches to use. The choice should give anintegration trajectory that follows significant energy in the scatteredwavefield and one that gives a well behaved approximation to thecontinuous integral. If the integral is performed as a weighted sum oversampled data, this means that the trajectory should be at leastpiecewise continuous.

One prior art technique as discussed briefly above is to use the firstarrival because it is easier to calculate and it is guaranteed to be acontinuous function. However, it has been pointed out that the firstarrival may contain very little energy and is thus not always a suitablechoice.

FIG. 3 shows a post-analysis seismic image in which the horizontal axisshows a distance (in feet) along the surface from the origin of thesurvey. The vertical axis shows the depth (in feet) beneath the surfaceof, in this case, a water layer that extends to 6000 feet. Thetraveltimes used in building this image were created using a finitedifference solution to the well known eikonal equation. The pre-stackmigration was performed using the first arrival at the seismic receiver.In this figure it can be noticed that many of the features are somewhathazy. In particular, there are two layers in the image locatedapproximately at 17,000 feet and 27,000 feet respectively. The clarityof these layers is very poor and there are significant parts of theimage of these layers where the layer is not evident at all. This isparticularly the case in the higher layer (17,000 feet) at distances ofapproximately 33,000 feet to 50,000 feet from the origin of the surveyand in the lower layer (at 27,000 feet) between 30,000 and 47,000 feetapproximately from the origin of the survey. There are two large saltbodies, the first located at between 6,000 and 14,000 feet in depth andbetween 10,000 and 37,000 feet from the origin of the survey. The secondsalt body is located between approximately 6,000 and 10,000 feet indepth and between approximately 45,000 and 75,000 feet from the originof the survey. As described above with reference to FIG. 2, the saltbodies can be responsible for refracted ray paths and this can result inlow energy levels in the first arrival. The images of the layers in thismodel are consequently poorly defined. In addition the grid-likefeatures at a distance between approximately 40,000 and 60,000 feet anda depth between approximately 9,000 and 17,000 feet are poorly defined.The analysis based on the first arrival also provides poor resolution ofthese comparatively small features.

Another prior art technique has been to choose the maximum amplitudearrival. This is the best single arrival to use as long as the estimatesof amplitude are correct. The reference cited previously showsimprovements in prestack depth imaging when maximum amplitudetraveltimes calculated by ray tracing are used.

FIG. 4 shows a seismic image covering the same distance and depth andderived from the same seismic model as the image shown in FIG. 3. Inthis case the maximum amplitude arrival was selected for pre-stackmigration and subsequent stacking. By comparison with FIG. 3 it can beseen that almost every feature in the image is sharper. The twohorizontal layers at approximately 17,000 feet and 27,000 feetrespectively have also been reconstructed rather more clearly than inthe image based on the first arrival. However, those parts of the twolayers identified in the discussion of the previous image are stillreconstructed somewhat unclearly. In particular, the lower layer (27,000feet) is reconstructed somewhat vaguely over the range betweenapproximately 32,000 to 44,000 feet from the origin of the survey. Inaddition, the small grid of features arranged at a distance betweenapproximately 40,000 and 60,000 feet are substantially clearer.

Unfortunately the amplitudes calculated in the asymptotic approximationare not necessarily a good approximation to amplitudes of the finitefrequency wavefield. The high frequency amplitudes are much moresensitive to the fine detail of the velocity model than the finitebandwidth wavefield. This has two results, first the choice of arrivalcan switch rapidly from one branch to another and secondly the choice ofbranch is very sensitive to changes in the model. Since we generallyonly have an approximate model, the choice made can be very different tothe choice that would be made in the true model. This clearlycompromises the accuracy of the analysis.

To ameliorate these shortcomings, the present invention provides a newcriterion for selecting the signal traveltime to be used for imaging.The traveltime that is selected is that associated with the ray that hasshortest physical length. This ray is often the highest amplitudearrival but this is not guaranteed. More importantly, it is almost nevera low amplitude first arrival associated with refracted energy. Inaddition, the physical length criterion is much more stable with respectto small fluctuations in the velocity model or to changes in the modelbecause it is an integral quantity as discussed above.

Unlike that for the first arrival, the traveltime from the shortest rayis not a continuous function. It is, however, piecewise continuous. Inaddition, because it is more stable, it tends to consist of a few largepieces with well defined boundaries between them. In contrast, theresult of using maximum amplitudes tends to be a much more “broken up”trajectory consisting of lots of small segments with rapid switchingback and forth between branches.

While a continuous function is desirable, the piecewise continuousfunction provided by the shortest ray has other advantages (e.g. it isassociated with more of the energy).

FIG. 5 shows a further image corresponding to those of FIG. 3 and FIG. 4and being derived from the same seismic model. In this instance,however, the pre-stack migration has been performed using the shortestray length rather than the ray corresponding with the first arrival orthe maximum amplitude arrival. By comparison with FIG. 3, it can readilybe seen that the image provided using shortest-ray prestack migrationresults in a much clearer image than selection of the first arrival. Theimprovement of this image over that derived using the maximum amplitudearrival is not quite so marked but an improvement nonetheless exists. Inparticular, the reconstruction of the layers at approximately 17,000feet and 27,000 feet respectively is improved. The reconstruction of theupper layer has improved clarity, especially between approximately32,000 feet and 50,000 feet from the origin of the survey. There is aneven clearer improvement over the image shown in FIG. 4 in the case ofthe lower layer (at 27,000 feet). The portion of this layer betweenapproximately 30,000 and 43,000 feet from the origin of the survey has adistinctly improved clarity. The continuity of both of these layersbelow the gap between the salt bodies is much improved. The advantagesin changing from maximum amplitude arrival to shortest path arrival isnot as great as that from the first arrival but it is generally clearthat the image everywhere is either just as good or it is improved.

In all of the migrated sections shown in FIGS. 3 through 5, the imagesof the curved events below 20,000 feet are artefacts due tomulti-reflection arrivals and should be ignored.

The selected raypath is used in prestack depth migration for two primarypurposes, to establish an acoustic pulse arrival time and to calculatean estimate of the amplitude of the arrival. These values may be used,for instance, to determine which input sample from a seismic traceacquired at the particular receiver location and associated with anacoustic pulse from the particular source location will be used tocalculate the output sample associated with the image point, as well asthe weight that will be applied to the input sample.

While FIGS. 3, 4, and 5 show two dimensional cross-sections of thesubsurface, they have been created using a three dimensional set ofseismic data. The prestack depth migration algorithm used to developthese sections is a true three dimensional method that does not assumethat the acoustic energy transmitted from the source to the receivertravels in the vertical plane connecting these points or that the imagepoint must reside in this vertical plane. This consideration of“out-of-plane” energy both increases the quality of the seismic imageobtained as well as increasing the complexity of properly processing thedata.

The selected raypath may also be used to update the seismic energypropagation model of the subsurface (i.e. the velocity model).

While the seismic data used to demonstrate the inventive method has beenconventional pressure-pressure (P-P) mode seismic data, the inventivemethod may be used in an identical manner with pressure-shear (P-S) modeseismic data (as well as shear-shear (S-S) mode and other seismic energytransmission modes) simply by providing a seismic energy propagationmodel of the subsurface that accounts for these alternative seismicenergy transmission modes.

The spirit and scope of the present invention is not restricted to thedescribed embodiments but encompasses any invention disclosed herein,explicitly or implicitly, and any generalization thereof.

1. A computer-implemented method of processing seismic data using aseismic energy propagation model of the subsurface, said methodcomprising the steps of: assigning seismic source, seismic receiver, andreflection point locations to said propagation model, identifying aplurality of alternative raypaths consistent with said propagation modelthat originate at said seismic source location, reflect at saidreflection point location, and terminate at said seismic receiverlocation, selecting a raypath from said plurality of alternativeraypaths having a shortest ray length, and utilizing said selectedraypath in subsequent seismic processing.
 2. A computer-implementedmethod as claimed in claim 1, in which said step of selecting a raypathcomprises calculating travel distances for each of said alternativeraypaths and comparing said calculated travel distances to determinewhich of said alternative raypaths has said shortest travel distance. 3.A computer-implemented method as claimed in claim 2, in which saidplurality of alternative raypaths are identified by a ray tracingprocedure.
 4. A computer-implemented method as claimed in claim 3, inwhich said propagation model is three dimensional and said plurality ofalternative raypaths are not coplanar.
 5. A computer-implemented methodas claimed in claim 4, in which said selected ray path is used tocalculate an acoustic pulse arrival time.
 6. A computer-implementedmethod as claimed in claim 5, in which said selected raypath is used toestimate an arrival amplitude.
 7. A computer-implemented method asclaimed in claim 6, in which said selected raypath is used to select asample from a seismic trace acquired at said seismic receiver locationassociated with an acoustic pulse from said seismic source location,said sample being used to create an image of said reflection point.
 8. Acomputer-implemented method as claimed in claim 7, in which saidselected raypath is used to calculate a weight to be applied to saidsample.
 9. A computer-implemented method as claimed in claim 8, in whichsaid propagation model models pressure-pressure (P-P) mode seismicenergy transmission in the geologic subsurface.
 10. Acomputer-implemented method as claimed in claim 8, in which saidpropagation model models pressure-shear (P-S) mode seismic energytransmission in the geologic subsurface.